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1. Introduction 

In recent years, ultracold atomic systems have served as a controlled and tunable 
toolbox for studying many-body quantum phenomena. The continuous tunability of 
the interaction by Feshbach resonances makes these systems ideal candidates to study 
the crossover from momentum space pairing in the Bardeen-Cooper-Schrieffer (BCS) 
theory to Bose-Einstein condensate (BEC) of fermions bound into molecules. This 
BCS-BEC crossover has been one of the most studied problems in recent experiments 
in both magnetic or optical traps m^lElEEllS and optical lattices [7]. 

Tuning across the Feshbach resonance, one traverses the whole range of the gas 
parameter kpOLs, where kp is the Fermi momentum and is the s-wave scattering 
length. The regime oi kF\as\ <C 1 (negative as) is described by Bardeen-Cooper- 
Schrieffer (BCS) theory. At kpas ^ 1 (positive as) fermions pair into bosonic 
molecules and form a Bose-Einstein condensate (BEC). At the microscopic scale, 
the BCS and BEC regimes are radically different; however, the macroscopic, and, 
in particular, critical behaviour is supposed to be qualitatively the same for the whole 
range of kpas'. the system undergoes the superfluid (SF) phase transition at a certain 
critical temperature. Separating the BCS and BEC extremes is the so-called unitarity 
point [kpas)^^ 0. It is worth noting that the unitarity regime is approximately 
realized in the inner crust of the neutron stars, where the neutron-neutron scattering 
length is nearly an order of magnitude larger than the mean interparticle separation 

m 

The Fermi gas at unitarity is a peculiar case of a strongly interacting system with 
no interaction-related energy scale: the divergent scattering length and any related 
energy scale drop out completely. This gives rise to universality of the dilute gas 
properties, in the sense that the only relevant energy scale left in the system is given 
by the density, n. Because of this universality one obtains a unified description of such 
diverse systems as cold atoms in magnetic or optical traps, Fermi-Hubbard model in 
optical lattices and inner crusts of neutron stars. 

The theoretical description of the Fermi gas in the BCS-BEC crossover regime is 
a major challenge, since the system features no small parameter on which one could 
build a theory in a rigorous way. The original analytical treatments were confined to 
zero temperature and were based on the extension of the BCS-type many-body wave 
function 9,. Most of the subsequent elaborations are also of mean-field type (with 
or without remedies for the effects of fluctuations) [El CH [El [El d HSl IIH ■ The 
accuracy and reliability of such approximations is questionable since they inevitably 
involve an uncontrollable approximation. 

Numerical investigations of the unitary Fermi gases are hampered by the sign 
problem, inherent to any Monte Carlo (MC) simulations of fermion systems |17l I18| . 
One way of avoiding the sign problem at the expense of a systematic error is the 
fixed-node Monte Carlo framework, which has been used to study the ground state 
The systematic error of the fixed- node Monte Carlo depends on the quality of 
the variational ansatz for the nodal structure of the many-body wave function and is 
not known precisely. Only in a few exceptional cases can the sign problem be avoided 
without incurring systematic errors. One of such cases is given by fermions with 
attractive contact interaction, for which a number of sign-problem-free schemes has 
been introduced [SOlEIlESEni- Fortunately, this system can be tuned to the unitarity 
regime. Still, despite a number of calculations at finite temperature |24l 1251 1^ . an 
accurate description of the finite-temperature properties of the unitary Fermi gas is 
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Figure 1. Diagrammatic series for the vertex insertion r(^, p) (heavy dot). 
Small dots represent the bare Hubbard interaction U, and lines are the single- 
particle propagators. 

missing. 

In the present paper, we simulate the Fermi-Hubbard model in the unitary regime 
by means of a determinant diagrammatic Monte Carlo method. By studying the dilute 
limit of the model, we extract properties of the homogeneous continuum Fermi gas. 
A brief summary of the main results has been given in Ref. ^7^. Here we provide a 
detailed description of the Monte Carlo scheme and methods of data analysis. We also 
report new results relevant to experimental realizations of the Fermi-Hubbard model 
in optical lattices and trapped Fermi gases. 

The Fermi-Hubbard model is defined by the Hamiltonian 

H = Ho + Hi, (1.1a) 
^^0 = ^(ek -M)Ck<TCk<T, (1-1&) 

kcr 

Hi = U^n^^n^l. (1.1c) 

X 

Here cj^^ is a fermion creation operator, Uxo- — cl^^C:^^, o" =T,i is the spin index, x 
enumerates sites of the three-dimensional (3D) simple cubic lattice with periodic 
boundary conditions, the quasimomentum k spans the corresponding Brillouin zone, 
ek = 2<^^^-^(l — cosfcao) is the single-particle tight-binding spectrum, a and t are 
the lattice spacing and the hopping amplitude, respectively, /z stands for the chemical 
potential and J7 < is the on-site attraction. Without loss of generality we henceforth 
set a and t equal to unity; the effective mass at the bottom of the band is then m = 1/2. 

In Sec. |2| we will study the two-body problem at zero temperature, show how 
the Hubbard model can be used to study the continuum unitary gas and investigate 
the functional structure of lattice corrections to the continuum behaviour. In Sec. 
Owe discuss the finite-temperature diagrammatic expansion for the Hubbard model 
(Sec. I3.1|l . and give a qualitative description of the Monte Carlo procedure to sum 
the diagrammatic series fSec l3.2|l . with details of the updating procedures given in 
Appendix [Appendix A| In order to extract the thermodynamic limit properties from 
MC data, we use finite-size scaling analysis described in Sec. 01 Section |31 gives an 
overview of the scaling functions describing thermodynamics of the unitary gas, and 
results are presented and discussed in Sec. El 

2. Two-body problem 

Consider a quantum mechanical problem of two fermions at zero temperature 
described by the Hamiltonian H1.1q|) - ()1.1c|I with /z = 0. The most straightforward 
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way to tackle this problem is within the diagrammatic technique in the momentum- 
frequency representation [^I^HI, which, in the present case, is built on four-point 
vertices, t/, with two incoming (spin-l and spin-|) and two outgoing (spin-l and spin- 
l) ends, connected by single-particle propagators. The scattering of two particles is 
then described by a series of ladder diagrams |29[ §16] shown in figure ^ Ladder 
diagrams can be summed by introducing the vertex insertion r(^, p), which depends 
on frequency ^ and momentum p. Since r(0, 0) is proportional to the scattering 
amplitude, the unitarity limit corresponds to r(^ — > 0, p ^ 0) — > oo. The summation 
depicted in figure pleads to 

T-\tp)^U-'+m,P) , (2.1) 
where n(^, p) is the polarization operator (the integration is over the Brillouin zone): 

m,P)- f T^TT ■ (2-2) 

It immediately follows from Eq. H2.1|l and (|2.2|l that the unitary limit corresponds 
to U = C/,, where 



U,^^-niO.O)^-l^l-. ,2.3) 



A straightforward numeric integration yields C/* ~ — 7.915i. 

In the limit of vanishing filling factor ^ for the many-body problem of Hl.la|l - 
1)1. lc|) the typical values of ^ and p are related to the Fermi energy ^ ^ ep v^^^ 
and Fermi momentum p ^ kp ^ v^l^ and are small compared to the bandwidth and 
reciprocal lattice vector, respectively. In zero-th order with respect to the lattice 
system is identical to the continuum one. Indeed, by combining H2.1|l . H2.2|l and (|2.3|l . 
we get 



r-^(e,p)= / 



dk 



1 1 



+ ep/2+k + ep/2-k 2ek_ 

and observe that for small ^ and p one can replace e(k) with fc^ and extend integration 
over dk to the whole momentum space with the result 



rc-o^.(e,P) = -^V^ + |^. (2.5) 

With this form of Fcont, and particle propagators based on the parabolic dispersion 
relation one recovers the continuum limit behaviour. 

Now we are in a position to treat the lattice corrections. The first correction 
should come from F, not from propagators, since only in F large momenta play a 
special role due to resonance in the two-particle channel. In the lowest non-vanishing 
order in ^ and p we have (the summation over repeating subscripts is implied) 

r-i^_ f g + (l/4)(9^£k/9fc,9fc,)p,p, 

^ 4(2^)3 44 ■ ('-'^ 

with the difference between the lattice and continuous model given by 
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J (27r)3 (fc72m)2 J^^ [2^f ^ - ( ■ ) 

[ dk 1/m /• dk i&'ejdk^dk^) 

J (2^)3 (fc72™)2 (2^)3 4 ■ I ■ J 

In the Umit of ^ ^ and p 0, we have T^^ « T^Tq^jj ^ kp ^ v^^^ , and 
— r^^j ^ k\ ^ v^l^ . Hence, the leading lattice correction is of the form 

r(e,P) ^ • ^ ^ ^ 

Incidentally, Eqs. H2.8|l and H2.9() hint into an intriguing possibility of completely 
suppressing the leading-order lattice correction by tuning the single-particle spectrum 
ek so that A = i? = 0. We did not explore this possibility in the present study. 



3. Determinant Diagrammatic Monte Carlo 

The diagrammatic technique employed in the previous section is not particularly 
convenient for numerical studies. In this section, we briefly review the Matsubara 
technique and then present a Monte Carlo scheme of summing the resultant 
diagrammatic series. 



3.1. Rubtsov's representation 

To construct a diagrammatic expansion for the model Hl.lfl|l - Hl.lc() we follow 
Refs. 1^ inni and consider the statistical operator in the coordinate — imaginary time 
representation. In the interaction picture we get: 



exp(-/3H) = exp{-PHo)Trexp(^- J\tHi{t)^ , 



(3.1) 



where (3 is an inverse temperature, Hi{t) — e'^^° Hie^'^^" , and %- stands for the 
imaginary time ordering. 

Expanding Eq. H3.1|l in powers of Hi , one obtains for the partition function: 



oc 



n— xi...x^ 



0<ri<r2<...<,3 : 



e ]^c|(xjTj)ct(xjTj)ct(xjrj)c|(xjrj) 



.(3.2) 



Expansion (|3.2|l generates the standard set of Feynman diagrams. Graphically, 
the diagrams are similar to those of Sec. |21 and consist of the four-point vertices, 
[/, connected by the single-particle propagators G^a\'^i — Xj,Ti — Tj]^,(3) = 
— Tr [Tre-~^^°c\{^iTi)ca{'y^jTj)\. The p-ih order of the perturbation theory is then 
graphically given by a set of (p!)^ possible interconnections of vertices by propagators, 
see figure 121 

The diagrammatic expansion H3.2() is unsuitable for the direct Monte Carlo 
simulation since it has a sign problem built in: different terms in the series have 
different signs — a closed fermion loop brings in an extra minus sign [2H1- The trick 
is to consider all diagrams of a given order p with the fixed vertex configuration 

5p {(xj,Tj), j = (3.3) 
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Figure 2. Diagrammatic series for the partition function. Upper line is the 
graphical representation of the series 13. 2i , lower line depicts Eq. 13.41 . Diagram 
signs arc shown explicitly. 

as one. This implies summation over the (p!)^ ways of connecting vertices by 
propagators. Upon summation, Eq. H3.2|l takes on the form |22|: 

oo 

Z = ^(-[/)P^detAT(5p)detA^(5p) , (3.4) 

p=0 Sp 

where 

E / ftdr,' (3.5) 

and A'^{Sp) are the p x p matrices buih on the single-particle propagators: 

Af^.(5p) = Gt'')(x,-x„T,-Tj), i,j = l,...,p. (3.6) 

For equal number of spin-up and spin-down particles detA^detA-'- = jdetAp, 
and the sign problem is absent. | Graphically, Feynman diagrams in this 
representation are just collections of vertices, see figure |5] For future use, we define 
the set of all possible vertex configurations H3.3(l by S*^^\ i.e., G^^^^ = {p, {Sp}}. 

The following two-point pair correlation function will prove useful: 

G2(xr;xV) = ( r.P(x, r)Ft(x', r') ) ^ 52(x^^ ^ (3^^^ 

<?2(xT;x'r') = Trr,P(x,r)pt(x',r')e-'3^, (3.8) 

where P(x, t) and pt(x',T') are the pair annihilation and creation operators in the 
Heisenberg picture, respectively: P(x, r) — Cx|(t)cxx(t). The non-zero asymptotic 
value of G2(xr;x'T') as |x — x'| — > 00 is proportional to the condensate density. 

Feynman diagrams for (72(xT;x'r') are similar to those for Z, but contain two 
extra elements: a pair of two-point vertices with two incoming (outgoing) ends which 
represent P(x, t) ( pt(x',r') ), see figure 13 The vertex configurations for the 
correlation function H3.7() slightly differ from those for the partition function l|3.3|l 
by the presence of the two extra elements: the configuration space for Eq. (|3.7() is 
6(G) ={p,{5p}}, with 

5p = {P(x,T),pt(x',r'), (x,,T,), j = l,...,p}. (3.9) 

J At half filling, the sign of U changes if the hole representation is used for one of the spin components. 
Hence, this method is also applicable to the half-filled repulsive Hubbard model. 



The Fermi- Hubbard model at unitarity 



7 



o + {T> + oo + ••• 
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Figure 3. Diagrammatic series for the correlation function 113.71 . Diamonds 
represent the two-point creation/annihilation operators P and Pt. 

The partially summed diagrammatic expansion for g2{ii.T]ii! t') is similar to Eq. (|3.4|) : 

oo 

g2(xT;x'T') = ^(-C/f ^detAf(4)detA^(4) , (3.10) 

where A'^{Sp) is a (p + 1) x (p + 1) matrix which differs from Eq. H3.6(l only by that it 
has an extra row iq and an extra column jo such that — G^a ^ (x^ — x, — r) and 

<, = G(°)(x'-x„r'-r,). 

Below we deal only with equal number of spin- up and spin-down fermions and (for 
the sake of brevity) suppress the spin indices wherever possible. We also peruse the 
generic notation (with superscripts) 'D{Sp) for p-th order terms of the diagrammatic 
expansions similar to ig^, e.g., V^^\Sp) = (-t7)2'|det A(5p)| . To simplify the 
notation we also omit superscripts if this does not lead to ambiguity. 

3.2. Diagrammatic Monte Carlo and Worm algorithm 

Equations (|3.4|) and H3.10|l have similar general structure of a series of integrals 
and sums with ever increasing number of integration variables and summations. In 
Refs. [21 122] it has been shown how to arrange a numerical procedure that sums 
such convergent series. To this end one considers the space of all possible vertex 
configurations 6 ( for the series ()3.4|l 6 = 6'^'^\ while for the series (|3.10|l 6 = 6'*^-' ), 
with the "weight" T>(Sp) associated with each element of the space. One then uses the 
Metropolis principle jSH] to arrange a stochastic Markov process which sequentially 
generates vertex configurations Sp according to their weights 'D{Sp), thus sampling 
the space S. In the course of sampling, one also collects statistics for observables in 
the form of MC estimators, see Sec. 13.31 

The stochastic process consists of elementary MC updates performed on vertex 
configurations Sp. The set of updates is problem-specific being restricted only by the 
requirements of (i) the ergodicity, i.e., given a particular diagram Sp it takes a finite 
number of steps to convert it into any other diagram S'^, and (ii) detailed balance, i.e., 
relative contributions of diagrams Sp and S'^ to the statistics is given by the ratio of 
their weights, VlSp) /ViS'^). The set of updates satisfying these requirements is not 
unique, the freedom is used to maximize the efficiency of simulations, as explained in 
detail in [Appendix A| 

In view of close similarity between the diagrammatic expansions (|3.4|l and (|3.10|l 
it is advantageous to construct a Monte Carlo process which samples these two scries 
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in a single simulation. This way one has access to both diagonal, e.g., energy, and ofF- 
diagonal properties, e.g., the superfluid response. An efficient way of performing such 
concurrent simulation is provided by the worm algorithm, which was originally devised 
for the worldline Monte Carlo simulations In the context of the diagrammatic 

determinant Monte Carlo, the generic worm algorithm principles imply the following. 
First, one works in the joint configuration space ©^'^^ 06^*^-* , accommodating diagrams 
l|3.3(l and H3.9|l . Second, all the updates are performed exclusively in terms of the two- 
point vertices P(x, r) and P^(x, r) — through their creation/annihilation, motion, and 
"interactions" with adjacent vertices. The worm-type updating procedures are further 
detailed in [Appendix A| Within the worm-algorithm framework, the configuration 
spaces 6^^^ and &^^^ are disjoint subsets of one extended configuration space. In 
what follows we will refer to them as Z-{ot "diagonal") and G-(or "off-diagonal") 
sectors of the configuration space. 

Formally, the extended configuration space corresponds to the generalized 
partition function 

- Z + CY, / dr / dr'g2(xr;xV) , (3.11) 

where the value of C is arbitrary — it controls the relative statistics of Z- and G-sectors 
and the efficiency of simulation. 



3.3. Monte Carlo estimators 

Suppose we have an observable X{q), which depends on a set of variables a, e.g. 
temperature and chemical potential. A MC estimator for the observable X is an 
expression which, upon averaging over the sequence of MC configurations converges 
to the expectation value of X{a). 

In accordance with Eq. (|3.11|) , the simplest worm-algorithm MC estimators are 

f 1, 5p e 6(^) , 



and 



Their MC averages are 



and 



^^,Jq Jo 



(3.15) 



where (. . .)mc denotes averaging over the set of stochastically generated configura- 
tions. In particular, for our purposes it will be quite useful that 



(,5(G)) 



MC 



CV / dr / dr'G2(xT;x'T') 
^^,Jo Jo 



MC 



(3.16) 
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The general rules for constructing an estimator for a quantity ^(o;) specified by 
the diagrammatic expansion 

X(a) = ^I?W(a;5p) , (3.17) 

are standard. We adopt a convenient convention: If the actual summation in 
H3.17|l involves only a subset So of vertex configurations — a typical example is an 
expansion defined within the Z-sector only — then we extend summation over the entire 
configuration space by simply defining V'^^^Sp ^ &o) = 0. If vertex configurations Sp 
are sampled from the probability density V^^^'^Sp) which comes from the expansion 
for the generalized partition function: 

ZH.(a) =^I?(^'*)(a;5p) , (3.18) 
then the MC estimator for X{a) is derived from 

(^) - ' (3-19) 

\" /MC 



as 



In what follows, by the estimator for a quantity x{a) — X(a)/Z{a) we basically mean 
corresponding function Q^^\a;Sp). 



3.3.1. Estimators for number density and kinetic energy We start with the estimator 
for the number density. The expectation value of the number density reads 

(3.21) 



2Tr4^(r)cxcr(r)e 



Z 

Here (x,r) is an arbitrary space-time point (the system is space/time translational 
invariant) and a is one of the two spin projections; the factor of 2 comes from the spin 
summation. The diagrammatic expansion of the numerator is similar to that for Z, 
with the diagram weight given by 

pt^H^p) = 2(-[/)P detB^+i(5p,x,r) detA;"(5j,) , (3.22) 

Here Sp G Ap'^{Sp) is a p x p matrix (|3.fci|) . and Bp_|_;^(5p, x, t) is a similar 

(p + 1) X {p + I) matrix with an extra row and a column, corresponding to the extra 
creation and annihilation operators in the numerator of l|3.21ll . respectively. This 
immediately leads to the following estimator H3.20|l for the number density 

We utilize the freedom of choosing (cr, x, t) to suppress autocorrelations in 
measurements. The density measurement starts with randomly generated (cr, x, r). 

The estimator for kinetic energy is derived similarly. One employs the coordinate- 
space expression for the kinetic energy in terms of hopping operators and deals with 
a slightly generalized version of Eq. (|3.21() : 
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The rest is identical to the previous discussion up to replacement Bp^j^(iSp, x, t) — > 
Bp_|_j^(iSp, X, g, r), since now the spatial position of the creation operator is shifted 
from that of the annihilation operator by the vector g. In our case, only the nearest- 
neighbor correlator (|3.24|) has to be computed. 

3.3.2. Estimator for the interaction energy The estimator for the interaction energy 

Tr ffi e~^^ 

{H,) = (3.25) 

is readily constructed using a generic trick of finding the expectation value of operator 
in terms of which the perturbative expansion is performed. Consider the Hamiltonian 
H{X) = Hq + XHi and observe that 

Tri7,e-^- = ---Tre-^-.---. (3.26) 

The differentiation of Eq. H3.2I) for Z — Z{X) and letting A = 1 afterwards is 
straightforward since the diagram of order p is proportional to A^. Hence 

Q^"'HSp) = -p-^pS^^'HSp) . (3.27) 

3.3.3. Estimator for the integrated correlation function Following the general 
treatment of Ref. [32], one can construct an estimator for the correlation function 
H3.7|l . In this work we just need to sum and integrate this correlator over all its 
variables (see Sec.^: 

K{L,T) = i(3L^y^y] / dr / dr'G2(x - x', r - r') , (3.28) 

and the estimator for K(L, T) is particularly simple: 

Q(^)(5p) = (/3L'*)-2c-i(5'°)(5p) . (3.29) 



4. Extrapolation towards macroscopic continuum system 

The MC setup discussed in SecOlworks in the grand canonical ensemble with external 
parameters {L, T, /i}. In order to extract the critical temperature of a continuum gas 
from MC data, one has to perform the two-step extrapolation, (i) Upon taking the 
limit oi L —^ oo one obtains Tdp), the critical temperature of a lattice system at a 
given chemical potential, and translates it into Tc{v) by extrapolating the measured 
filling factor to the infinite system size: v = T = T^^), L — > oo). (ii) The 

extrapolation towards the continuum limit is then done by taking the limit of ^ 0. 
The latter procedure is based on results presented in Sec. [3 

The finite-size extrapolation is performed by considering a series of system sizes 
Li < L2 < . . . . At the critical point the correlation function (|3.7() decays at large 
distances as a power-law: 6*2 (x — x', t — t') cx |x — x'|^(^+''\ where 77 is the anomalous 
dimension '35'. Since we expect the transition to belong to U(l) universality class, we 
take rj « 0.038. If one rescales the summed correlator H3.28II according to 

R{L,T) ^ L'^+'^K{L,T) , (4.1) 

the corresponding quantity is supposed to become size-independent at the critical 
point, i.e. the crossing of R{Li,T) and R{Lj,T) curves can be used to obtain an 
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estimate Tj^^^l for the critical temperature T^ip) Indeed, for temperatures in 

the vicinity of the critical point the correlation length diverges as ^corr oc 
where t = {Tdp) — T)/Tc{fJ,), and « 0.671 for the U(l) universality class. In 
the renormalization group (RG) framework , the finite-size scaling of the rescaled 
correlator R obeys the relation 

i?-/(a;)(l + cL- + ...) , (4.2) 

where x = {L/^corvY^"'' is the dimensionless scaling variable, f{x) is the universal 
scaling function analytic at x = 0, c is a non-universal constant, oj « 0.8 is the 
critical exponent of the leading irrelevant field [37], and dots represent higher-order 
corrections to scaling. If the irrelevant field corrections were not present, all R{Li, T) 
curves would intersect at a unique point, Tdp). Expanding Eq. H4.2|) up to terms 
linear in t one obtains for the crossing T^.^l. 

_ const {Lj/UY - 1 

To employ Eq. (|4.3|l one performs a linear fit of the sequence of Tl^^l against the 
right hand side of Eq. (|4.3|l for several pairs of system sizes. The intercept of the 
best-fit line yields the thermodynamic limit critical temperature Tc{fi). Note that if 
the universality class is not U(l) and the values of rj, v^, to are different, or system 
sizes are too small to justify the scaling limit, the whole procedure fails. Hence, the 
adopted scheme of pinpointing Tc features a built-in consistency check. 



5. Thermodynamic scaling functions of a unitary Fermi gas 

As has been noted in Sec. ^ the only relevant microscopic energy scale in the 
continuum unitary fermi gas is given by the Fermi energy ep = nh^n^/'^ /m, where 
K = (37r^)^/^/2 for a two-component Fermi gas. Therefore, as it was first noticed in 
Ref. j38| . all thermodynamic potentials feature self-similarity properties and can be 
written in terms of dimensionless scaling functions of the dimensionless temperature 
X = T /ep. All scaling functions are mutually related; it is sufficient to know just one of 
them to unambiguously restore the rest. Apart from the shape of scaling functions, the 
self-similarity at the unitary point is identical to that of a non-interacting gas Fermi 
gas |39| . including functional relations between different thermodynamic potentials. A 
derivation of scaling functions and relations between them can be found in Ref. |38| . 
In this section, we render the scaling analysis in a form convenient for our MC study. 

In terms of the dimensionless chemical potential y = /i/sF, the dimensionless 
equation of state reads y = f^{x). The function can be calculated numerically. 
Another quantity which is also available in our simulation is the dimensionless energy 
per particle E/{NeF) — fE{x). The scaling relations for other thermodynamic 
quantities are defined likewise. For instance, the entropy and pressure read S/N = 
fs{x), and P/{n£p) = fp{x). Though fg and fp are not directly calculated in our 
simulation, and we can relate them to fp- It is also important to relate to fp, 
since this yields a consistency check for the numerical results. 

To establish desired relations, we start with the scaling relation for the Hclmholtz 
free energy F/{Nep) = fp{x) which in canonical variables reads 

F(r, N,V) =7 fp (t/7(A^/F)2/3) {N/Vf/^ N , (5.1) 
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where 7 = nh? /m. Then, for the entropy and pressure we have (the prime stands for 
the derivative) 

fs = -f'p, (5.2) 
fp = {2/i){fF- f'px) . (5.3) 

The expression for energy in terms of fp is 

fE = fF-f'pX. (5.4) 

We thus see that 

fp EE (2/3)/£ . (5.5) 

On may also consider Eq. H5.4|l as a differential equation: 

fF-fpX^fE (5.6) 

to be integrated with respect to fp from x — 00 down to finite x, taking advantage of 
known asymptotic behaviour of fp and fp for the weakly interacting two-component 
gas. 

Now, we note that from the general thermodynamic relation E — —PV+TS+fj,N 
it immediately follows that 

fE = -fp + xfs + , (5.7) 
which, in turn, leads to the following relations 

fs = mfE_A , (5.8) 

X 

fp = f^- {2/3) fp , (5.9) 

that allow one to express fs and fp functions through numerically available functions 
fp and /p. Another useful relation is 

f'^ = ^^"(^/^^^^ , (5.10) 

X 

which allows to extract fp directly from fp and J^, and thus provides a simple check 
for the data consistency: The result H5.9() for the fp curve should be consistent with 
the the derivative deduced from H5.10|l . 
By integrating equation H5.6|l we get 

fp{x) = Cqx - ^x\nx ~ X f ^ ) da;o ■ (5-11) 

^ J X K'^ Xq Xq / 

Here we took into account the asymptotic ideal-gas behaviour of fp: 

3 

fp{x) ^2"'' a; ^ 00 , (5.12) 

and introduced the corresponding term into the integral to render the latter 
convergent. The free constant of integration, Co, can be restored from H5.9|) - (|5.12|l 
combined with the asymptotic ideal-gas behaviour of 

fti{x) — > — -a:;ln( — a;) at a; — ^ 00 . (5.13) 

2 V ZTT / 



The result is 



Co = - 1. (5.14) 



K I 
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Note that if higher-order terms in the asymptotic {x — s- oo) behaviour of fsi^) are 
also known, then H5.11|) can be used to estabhsh the corresponding corrections for 
fpix) and other scahng functions. For instance, it has been found in Ref. that, 
as a; ^ oo 

/^.(-)-l--l(-)'^'^- (5.15) 
2 8 Vk/ y/X 

In accordance with 15.11|l . this imphes (x — )■ oo) 

3 3 /7r\3/2 1 

frix) Cqx- -x\nx- - [-) (5.16) 
2 4 Vk/ yjx 

, , ^ 3 / K \ 3 /7r\3/2 1 

FinaUy, the scahng functions Wo and Qa defined in Ref. [HH! are related to Je and 
as foUows: 

40 jEix) 
a;5/2 ' 



Wo(^(x)/2:) EE TT^^, (5.18) 



^0 (-//mW) - ^ ■ (5.19) 

5.1. Trapped Fermi gas 

So far we have considered the uniform Fermi gas, while in experimental realizations 
m m El 01 |S1 ini one has to deal with the parabolic trapping potential. The standard 
procedure, especially in systems with short "healing length" is to use the local density 
approximation (LDA), i.e. to replace the chemical potential with its coordinate- 
dependent counterpart /i(r) — fJ- — V{r). This procedure can be easily combined with 
MC results as follows. We introduce the dimensionless variable u = fi/T ~ f^{x)/x 
and define the scaling function for the number density as w„ = x^^/"^ . This is 
equivalent to the parametric {u{x),Wn{x)} dependence of n on u. The scaling 
functions for other thermodynamic quantities are defined in a similar manner, e.g. 
we{u) = fE{x{u)) for the energy, and ws = fs{x{u)) for the entropy. Within LDA 
u acquires the coordinate dependence u{r) = (p — V(r))/T which translates into the 
density profile 7i(r;/i, T) = Wn{u{v))[mT / nh?Y/'^ . Likewise, other thermodynamic 
functions are to be understood as local, coordinate-dependent quantities. 

Consider the case of N particles in a cigar-shaped parabolic trap, characterized 
by the axial and radial frequencies uj\\ and The characteristic energy in this case 
is Ep — (3iV)^/'^ft([x)j|[x)j^)^/'^, which would coincide with the Fermi energy for a non- 
interacting gas in the trap. Note, that we denote the Fermi energy in the trap by an 
upright capital Ep in order to avoid confusion with the uniform system Fermi energy 
£f- 

By integrating over the radial coordinates one obtains the axial density profile 
na{z) {z is an axial coordinate) in the form 

^ i2T/EEfl^ ( , \ 
N 7rL|| \T 2{T/EF)L\j' ^ ' 

where L|| — A^/^(3A^)^/^/|| , the aspect ratio A = a;x/tj||, the oscillator length 
Z| = h/mLu^\, and 

Wn{p) = / Wn{u)du. (5.21) 
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By integrating Eq. ()5.2U|I with respect to z, one finally relates chemical potential 
to temperature: 

where 

Wn{p)= / Wn{p-q'^)dq. (5.23) 

JQ 

Obtaining the temperature dependence of thermodynamic functions for a non- 
uniform system within LDA is also straightforward. For the total energy of a cloud 
Etot we obtain 

M = ? (|^)7;>pf (5.24) 

and likewise for the entropy: 



_ 16 



3 i-^/T 

dpi dqws(p- q^)wn(p~ q^)- (5.25) 



6. Results and Discussion 



We performed simulations outlined in previous Sections for filling factors ranging 
from 0.95 down to 0.06 with up to about 300 fermions on lattices with up to 16^ 
sites. The typical rank of determinants involved in computations of acceptance ratios 
(Sec. [Appendix A| | and estimators (Sec. 13. 3() is up to M ~ 5000. Since we only need 
ratios of determinants, we use fast-update formulas |22j to reduce the computational 
complexity of updates from down to AP . 

We validate our numerical procedure by comparing results against the exact 
diagonalization data for a 4 x 4 cluster and other simulations of the critical 

temperature at quarter filling v — 0.5 [421 143 | and v = 0.25 In all cases we find 
agreement within statistical errors of a few percent. 



6.1. Critical temperature 

Figure01shows a typical example of the finite-size analysis outlined in Sec. 01 Despite 
the fact that numerically accessible system sizes are quite small, figure 0] (and 
similar analysis for the whole range of filling factors) supports expectations that the 
universality class for the SF phase transition is U(l). The finite-size analysis allows 
us to pinpoint the phase transition temperature to within a few percent. 

Shown in figure El is the dependence of the critical temperature on the lattice 
filling factor. The critical temperature is measured in units of the Fermi energy, as is 
natural for the unitarity limit. We define the Fermi momentum for a lattice system 
with filling factor v as kp — (37r^i^)^/^ and the Fermi energy ep = kj^, as those of a 
continuum gas with the same effective mass and number density n = v. 

It is clearly seen that the presence of the lattice suppresses the critical temperature 
considerably, nearly by a factor of 4, depending on the filling factor. Strong dependence 
of Tc on V is in apparent contradiction with Ref. |24| . which claims weak or no in- 
dependence. This disagreement might be due to the difference in the single-particle 
spectra et used: Ref. j24) employs the parabolic spectrum with spherically symmetric 
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C F 




Figure 5. The scaling of the lattice critical temperature with filling factor 
(circles), u = 1 corresponds to the half filling. The errorbars are one standard 
deviation. The results of Ref. _42 43 at quarter filling and v = 0.25 are also 
shown for a comparison. See the text for discussion. 

cutoff, while we use the tight-binding dispersion law. Indeed, Eqs. (|2.8() - (|2.9() indicate 
that a particular choice of £k does influence lattice corrections to Tc, which may even 
have different signs for different £k- 
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Figure 6. The fit residues for the best-fit hne of figure O plotted versus h>^^^ 
(circles). Zero level is shown by the horizontal line, the blue dashed lines are 
linear fits to the data points for u < 0.5 and u < 0.35, respectively. 



It is also clear from figure |S1 that close to half- filling Tc is essentially constant, 
as expected (see, e.g. 01]). The predicted ^ v^^"^ scaling H2.HJ|I sets in at about 
u ^ 0.5. We thus use a linear fit Tc{v)/£p{i') — Tc/sp — const • v^^^ to eliminate 
lattice corrections in the final result. Such fitting procedure results in the best-fit line 
given by 0.152(7) — 0.13(2)j/^/'^. We further analyze the fit residues in order to estimate 
the effect of the sub-leading lattice corrections which are expected to be proportional 
to iP'/^ . As shown in figure|^ such corrections, if any, are smaller than the uncertainty 
of the u^/^ fit. 

This analysis yields the final result Tc/ep — 0.152(7) for the continuum uniform 
gas, which is noticeably below the transition temperature in the BEC limit Tbec = 
0.218eF- Various approximate analytical treatments led in the past to either above 

[iniiTiiniiTni, or beiow [niiniiini ?bec- 

Is is instructive to compare our results for Tc to other numerical calculations 
available from the literature. The simulations of Ref. yield Tc = O.Qbep, but at 
the value of the scattering length which has not been determined precisely. This result 
most probably corresponds to a deep BCS regime, where the transition temperature 
is exponentially suppressed. Lee and Schafer report an upper limit Tc < O.lAep, 
based on a study of the caloric curve of a unitary Fermi gas down to T /ep — 0.14 
for filling factors down to v = 0.5. The caloric curve of Ref. [221 shows no signs of 
divergent heat capacity which would signal the phase transition. This upper limit is 
consistent with Tc(t^ — Q.b)/£p ~ 0.054, see figureEl 

The Seattle group has performed simulations of the caloric curve and condensate 
fraction, no, of the unitary gas, Ref. Using "visual inspection" of the caloric 

curve shape the critical temperature was estimated in Ref. [22 to be Tc ~ 0.22(3)eF- 
Unfortunately, the authors did not perform the finite-size analysis and v ^ 
extrapolation. The overall shape of the caloric curve seem to be little affected by the 
finite volume of the system. This is hardly surprising since even in the thermodynamic 
limit E{T) and its derivative di?/dT are continuous at the transition point. These 
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0.05 0.10 0.15 0.20 0.25 0.30 7/5^ 

Figure 7. The finite-size scaUng of tlie condensate fraction data from Ref. 1241 . 
Raw data points are rescaled similar to Eq. 14. li by the L^^'^ factor. Shaded 
vertical strips represent results for Tc/sp of this work and Ref. 1241 . respectively, 
solid lines are drawn to guide an eye. 

properties also make it hard to use non-quantitative measures for reliable estimates of 
critical parameters from the E(T) curve. On the other hand, the condensate fraction 
which has singular properties at Tc does show sizable finite-size corrections, see figure 1 
of Ref. [3] . At this point we note that scaling of the condensate fraction is identical 
to that for K[L, T). In figure [3 we plot the data of Seattle's group as uqE^^^ versus 
temperature. The intersection of scaled curves turns out to be inconsistent with the 
estimate for Tc derived from the caloric curve inspection. 

6.2. Thermodynamic functions 

The filling factor dependence of thermodynamic quantities is similar to that of T^. 
Figure |H1 displays the behaviour of energy and chemical potential along the critical 
line T = Tdv). The extrapolation towards v yields for the continuum gas 

E/iNep) = 0.31(1) (r = T,) , (6.1) 

^l/eF = 0.493(14) (T = . (6.2) 

The numerical values for other thermodynamic functions at criticality can be easily 
restored using the formulas of Sec. \E\ 

In order to elucidate the thermodynamic behaviour of the unitary gas, we 
performed simulations for a range of temperatures T > T^. Shown in figure El are 
the simulation results for energy and chemical potential for the continuum gas as 
functions of temperature. Each point was obtained using data analysis similar to that 
depicted in figure|Sl In the high-temperature region we simulated up to 80 fermions on 
lattices with up to 32^ sites. In this region, the condition v \s necessary but not 
sufficient for extrapolation to the continuum limit, for it is crucial to keep temperature 
much smaller than the bandwidth: T <C 6i. 
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Figure 8. Energy (loft-hand panel) and chemical potential (right-hand panel) 
dependence on the filling factor along the critical isotherm T = Tc(v). Dots are 
the MC results, dashed lines are the linear fits. 



As can be seen from figure |H1 our results for both energy and chemical potential 
approach the virial expansion 40, as T/ep oo. For T/ep ^ 0.5 our data are not 
far from the curve of Ref. [23]. Though we do not have data points for T < Tc there 
is still a reasonable agreement even at Tc with the T ^ fixed-node MC values |19j . 
In this region, our results are consistent with a very weak dependence of energy and 
chemical potential on temperature, and the numeric values of both are consistent with 
the experimental results O 01 |S] . 

Using Eq. (|5.9I) and data from figure we deduce the dependence of free energy 
on temperature, see figure ITUl We also use Eq. H5.10|l to make sure that our MC data 
for energy and chemical potential are consistent with each other. 

6.3. Trapped gas 

As discussed in Sec. l5.1l the thermodynamic functions for the uniform case can be used 
for analysis of experimental system within the local density approximation. In this 
section we report preliminary results of our ongoing study of trapped gas experiments. 

It starts with the interpolation procedure which produces continuous functional 
behaviour for thermodynamic functions, consistent with the discrete set of simulated 
points. We use a piecewise-cubic ansatz with a smooth crossover to the virial 
expansion, Eq. (|5.1t)|) . for the free energy. Temperature dependence of both energy 
and entropy are then deduced using numerical integration of Eqs. H5.24|l and (|5.25|) . 
respectively. 

As the trapped gas is cooled down, the superfluidity first sets in at the centre of 
the trap, where the density is the highest. Equation (|5.22() can be used to pinpoint 
this onset temperature: at T = Tc, n/T = {iJ- /£%')/ {T / e^p) , where is the Fermi 
energy of the uniform gas with the density equal to the density at the trap centre. 
Using Tc/e%. = 0.152(7) and Eq. (|n3, one obtains Tc/Ep = 0.20(2). We quote here a 
conservative estimate for the uncertainty, which incorporates both the uncertainty of 
the critical temperature itself, and a systematic uncertainty which stems from restoring 
the continuous functional dependence of the chemical potential out of the finite set of 
the Monte Carlo calculated points with finite errorbars. 

Experimentally, the temperature of the strongly interacting Fermi gas is not easily 
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Figure 9. The temperature dependence of the energy per particle (upper 
panel) and chemical potential (lower panel) of the unitary Fermi gas. Red circles 
are the MC results, black dotted lines and blue dashed lines correspond to the 
Boltzmann and non- interacting Fermi gases, respectively, the dot-fashed lines are 
the asymptotic prediction of Ref. 1401 (plus the first virial Fermi correction), black 
triangles are the MC results of Ref. )24l . and the purple stars denote the ground- 
state fixed-node MC results 1191 . 



accessible. On the contrary, thermometry of the non-interacting Fermi gases is well 
established. In the adiabatic ramp experiments one starts from the non-interacting gas 
at some temperature [in units of Fermi energ y] (T/Tf) , and slowly ramps magnetic 
field towards the Feshbach resonance 0, thus adiabatically connecting the system 
at unitarity to a non-interacting one. Assuming the entropy conservation during the 
magnetic field ramp, equation (|5.25|l can be employed for the thermometry of the 
interacting gas: by matching the entropy of a non-interacting gas at the temperature 
(T/Tp)° with the entropy calculated via Eq. H5.25|l one relates the initial temperature 
(before the magnetic field ramp) to the final temperature (after the ramp). We find 
that the onset of the superfluidity corresponds to {T/Tp) = 0.12 ± 0.02 (again, we 
quote here the most conservative estimate for the errorbar). This value seems to be 
somewhat lower than the value suggested by the experimental results 0. Nevertheless, 
given the level of the noise in figure 4 of Ref. B the consistency is reasonable. 

An alternative thermometry can be build on recent advances in the experimental 
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Figure 10. Free energy versus temperature. Red dots are the MC data, and 
dashes represent the errorbar range for derivative of free energy, calculated via 
Eq. 15.101 . Black triangles are MC results of Ref. I24| . purple star denotes the 
ground-state fixed-node MC result 1191 . black dotted line shows the Boltzmann 
gas curve, and the blue dashed is the asymptotic prediction of Ref. 1401 . 



technique Pl E] which made it possible to directly image the in situ density profiles 
of the interacting system. Such density profiles can be directly fit to Eq. (|5.20|l . 
which gives the shape of the cloud depending on /i/T and T/Ep. By relating the 
chemical potential to the temperature using Eq. (|5.22|) . one is left with only one 
fitting parameter, T/Ep (apart from a trivial fitting parameter zq which accounts for 
the overall shift of the cloud image off the trap centre). 

As an illustrative example of such procedure we have analyzed the experimental 
density profiles measured by the Rice's group @], as depicted in figure ITTl From 
this analysis we deduce the upper bound for the temperature in the experiments 0] 
T < O.lEp, which is consistent with the results of the measurements of the condensate 
fraction 0. Since in the experiments [20 the gas is very degenerate, one is able to 
put an upper limit on temperature only. 

Note that if the temperature is known from, e.g. the adiabatic ramp experiments, 
Eqs. H5.20|l - H5.22|) must reproduce the cloud shape without free parameters (apart from 
zo)- 

7. Conclusions 

We have developed a worm-type scheme within the systematic-error-free determinant 
diagrammatic Monte Carlo approach for lattice fermions. We applied it to the 
Hubbard model with attractive interaction and equal number of spin-up and spin-down 
particles. At finite densities, the model describes ultracold atoms in optical lattice. In 
the limit of vanishing filling factor, — > 0, and fine-tuned (to the resonance in the two- 
particle s-wave channel) on-site attraction, a universal regime sets in, which is identical 
to the BCS-BEC crossover in the continuous space. In the present work, we confined 
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Figure 11. Axial density density profiles: experimental data (dots) is from the 
data contained in figure 3 of |^. The full red line is calculated via Eq. 15. 211 - 15.221 
with T/Ep = 0.03 (correspondingly, T/e^p^ = 0.02), dashed blue line corresponds 
to T/Ep = 0.16 (T/e^' = 0.1), and dot-dashed green line is for T/Ep = 0.22 

(r/e^°' = 0.16). In both cases we allowed for a horizontal displacement of a whole 
curve in the range of [2:0! < 20/im. See the text for discussion. 

ourselves to a special value of the on-site interaction, U = ~ — 7.915t, corresponding 
to the divergent s-scattering length. At [/ = [/» and 1/ ^ 0, the system reproduces the 
unitary point of the BCS-BEC crossover. The unitary regime is scale- invariant: all 
thermodynamic potentials are expressed in terms of dimensionless scaling functions 
of the dimensionless ratio T/ep (temperature in units of Fermi energy). We obtained 
these scaling functions by extrapolating results for the Hubbard model to — > 0. 
For the critical temperature of the superfluid-normal transition in the uniform case 
we found Tc/ep — 0.152(7). Our results form a basis for an unbiased thermometry of 
trapped fermionic gases in the unitary regime: In particular, we found (within the local 
density approximation) that for the parabolic confinement, the critical temperature in 
units of the characteristic trap energy Ep is T/Ep = 0.20(2). For the experimentally 
relevant case of an isentropic conversion of a gas from the non-interacting regime 
to the unitary regime we find that the onset of the superfluidity corresponds to the 
initial temperature (before the magnetic field ramp) (T/Tf)° = 0.12 ± 0.02, which 
is reasonably consistent with the experimental result (HI, to within the experimental 
noise. 
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Appendix A. Updating procedures 

The generic diagrammatic MC rules for constructing updates are as follows 1^ . Let an 
update B transform a diagram ^{Sp) into diagram Here configurations Sp and 

iS^ may have different order and may or may not feature the "worm" two-point vertices. 
The update B involves two steps. First, a modification of the configuration is proposed, 
with some probability density W{Sp — > S'^) for new continuous/discrete variables. 
There are no strict requirements fixing the form of W{Sp S'^), it should rather 
be chosen on physical grounds to maximize the efficiency of the algorithm, and be 
simple enough to allow analytical evaluation of the normalization integral. Then, the 
update is either accepted, with probability P^^^^/ , or rejected. The complimentary 
update C, which transforms T>{S'q) into 'D{Sq) proposes the modification with the 
probabiUty density W{Sp ^ S'^), which, in principle, may differ from W{Sp S'g), 
and is accepted with the probability Pg^^^s/^. For a pair of updates B and C to be 
balanced, the acceptance probabilities must obey the Metropolis relations 

^5,^5^ =min(l,7^), (A.l) 

P5,^5^ =min(l,l/7^) (A.2) 

where 7^ is a solution of the detailed-balance equation 

nw{Sp ^ s'^)v{Sp) = w{Sp ^ s'^)v{s'^). (a.3) 

If the probabilities of selecting updates B and C are not equal, they must also be 
included into the definition of W{Sp S'^) and W{Sp ^ Sg), respectively. 



Appendix A.l. The "Diagonal" updating scheme 



The minimal ergodic set of updates, as suggested in |521l3nij consists of just one pair of 
"diagonal" updates which increase/decrease the number of vertices in a configuration 
by one: in the forward update one adds a vertex at randomly selected point in the 
X (3 hypercube; in the reverse update one removes a random vertex from the 
configuration. In the simplest version, the probability density for selecting a point is 
uniform, i.e. one selects a particular lattice site with the probability p(xnow) = L^'^ ^ 
and a temporal position with the probability density w{t^^cvj) = ■ A vertex to be 
removed is selected at random out of p vertices present in the configuration. Hence 
the Metropolis acceptance ratio function IjA.Sp is given by 

n CT, /M 1-)., I 1 1 ^ 

T^add = 



det A(5j,+i 



det A{Sp 



{-U)I3L'^ 



P 



1 



(A.4) 

Here 5p+i = (xnowj '''new) U Sp is the configuration with an extra vertex (xngw; ''new)- 
The acceptance ratio to remove a vertex is 

det A(5p_i^ 



7^r 



det A{Sp 



{-U)l3Ld 



(A.5) 
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Figure Al. Multi-ladder diagrams in deep BEC regime (A), and close to 
unitarity (B). 



Appendix A. 2. Worm-type updating scheme 

The diagonal updating scheme described in the previous section is highly inefhcient 
in the dilute regime — the regime we are mostly interested in. To assess the efficiency 
of the "diagonal" scheme consider first a low-density gas in the deep BEC regime. 
In this limit, all fermions are paired into composite bosons (apart from exponentially 
rare fluctuations). Inter-boson interactions occur only via rare collisions. In the 
diagrammatic language, the existence of compact composite bosons is reflected by the 
fact that the dominant diagrams are the multi-ladder ones (see figure IXTK '). with the 
typical x-span of each ladder being of the order of L, and the typical r-span of the 
order of (3. Each ladder of such a diagram in fact mimics a world line of a composite 
boson. 

The analogy between ladder parts of multi-ladder diagrams and world lines of 
composite bosons helps elucidating the structure of typical diagrams at resonance. 
Indeed, at the two-particle resonance composite bosons exist only virtually, but at the 
same time the energy cost of creating a boson is zero. Thus, multi-ladder diagrams 
still dominate, but the x- and r-span of individual ladders becomes finite and related 
to the inter-particle distance, see figure lAlB . Resonant ladders also become more 
"dilute" since the average distance between vertices diverges, but they still contain 
many vertices since the resonance formation involves distances much smaller than 

Proposing diagrams which disregard the nature of the resonant interaction results 
in small acceptance ratios. To significantly reduce this type of slowing down, we resort 
to a worm-type updating scheme which allows one to, loosely speaking, directly "draw" 
the multi-ladder diagrams. 

Appendix A. 2.1. Worm creation/annihilation updates These updates are the 
necessary ingredient of the worm-type scheme, since they switch between the diagonal 
H3.3() and off-diagonal H3.9|l diagrams. These updates transform diagrams as follows: 

6(^) 9 5p = {..., (xi,Ti), ...}t^ (A.6) 

{. . . , (xi, n), P(x, r), Pt (x', r'), . . .} = 5p e 6(«) . (A.7) 

To create a pair of two-point vertices P(x, r) and P^(x', t'), one first selects x and 
r uniformly over the spatial lattice and (0, /?) interval, respectively. One then selects 
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x' among the lattice sites in the spatial cube with the edge length I centered around x 
with equal probabilities, and selects r uniformly on the interval (r — Ar/2, t + At/2). 
The values of I and At are arbitrary. The overall probability density is thus given by: 

The reverse update attempts to remove P(x, t) and P(x', t') provided |xq, — xj^| < 1/2 
{a = 1, 2, 3) and |t — t'| < At/2, as prescribed by the balance requirement. 
The solution to the detailed balance equation ljA.31) then reads: 



n = 



det A(5p) 



det A{Sp) 



(A.9) 



In order to maximize the efficiency of the numerical scheme, it is advantageous 
to have the acceptance probabilities (jA.l|l - ljA.2|) to be of the order unity. On the 
other hand, the explicit macroscopically large factors in the right-hand side of Eq. 
IjA.Qj) render the acceptance probabilities macroscopic. The use of the extra weighting 
factor C is now clear: by choosing 1/C = f3 L'^ AtI'^ / ( one removes undesired factors 
from the right-hand side of Eq. I|A.9|) . The remaining freedom of choosing C has to be 
employed to further fine-tune the acceptance probabilities. The required values of C 
are typically of the order of unity. 



Appendix A. 2. 2. Adding/removing vertices within the worm framework These 
updates are central for the method since they change the diagram order. The main 
idea of these updates is: given a pair of "worm" two-point vertices P and P^, use, 
say, P^' (x, t) as a tip of a magic pen to add or remove vertices from the diagram (note 
that it suffices to use either P or P^ as a dynamical variable, the choice being only a 
matter of taste): 

5, = {..., (X1,T1), pt(x,T), ...}^ ^ (A.IO) 

{. . . , (Xi,Ti), (x,t), pt(Xnew,T-now), ■ • ■} = S^+l ■ (A.H) 

We have employed two versions of these updates. 



High-density version In the most simple yet useful version of a forward update, one 
selects Tnew and Xnow uniformly in the interval (t — At/2, t -|- At/2), and the spatial 
cube of the edge length I centered at x, correspondingly. In the reverse update one 
selects a vertex to be removed at random out of m choices, where m is the number of 
vertices (xi,Ti) such that \{xi)a — (xnew)a| < 1/2 and It— t„cw| < At/2. The solution 
to the detailed balance equation then reads 

2 



n 



det A{Sp) 



(-J7)AW« 



(A.12) 



Low-density version The strategy just described works well close to half-filling. In the 
low-density regime the multi-ladder diagrams dominate and the above scheme is not 
optimal for a number of reasons: (i) it does not respect the diagram structure which 
should feature (at least locally) well-defined vertex chains; (ii) the probability density 
according to which the value of Tnew is proposed should favour shifting P'l' towards 
smaller t-s. To see this, consider two particles in vacuum: the matrix elements of 
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P^(x, Ti)P(x, simply equal zero for T2 < ti.]; (iii) both large and small values of 
6t = T — Tnow have to be accounted for. Consider two fermions in vacuum again. The 
local probability density for a shift of length 5t and x^cw = x is proportional to the 
free diffusion propagator: w{St) oc (5t~^/^, which is singular as St 0. However, the 
mean value of St for this distribution, J w{St)5t d{5T), diverges at the upper limit. 
These properties of resonant ladders are ignored when shifts are proposed uniformly 
over some ±Ar interval. 

The requirements (ii) and (iii) are best met if, when adding a vertex, a new 
position of is proposed according to the probability density [cf. Eq. (|A.lip ] 



w(xncw-x; St) cx 



|G(")(x- 
0, 



if 



> T 



otherwise. 



(A.13) 



Since the free single-particle propagator is known only numerically, we use Eq. PTT3)l 
on a uniform mesh with some step a: St-j = aj. That is, we tabulate the values 
Wj = w(y; Stj) prior to the start of computations; we then select y and Stj according 
to the (normalized) distribution w{y;STj). The proposed new position of P^ is then 
chosen as Xnow = x + y, and Tncw = t + Stj + A, where A is uniform over the interval 
[—a/2, a/2]. We find that a ^ 1/5U produces good results. 

In order to meet the requirement (i), the notion of the closest — in terms of a 



certain distance — neighboring vertex is introduced: the reverse update S' 



P+i 



deterministically deletes the closest neighbor of P^ . The detailed balance than requires 
that the forward update should be automatically rejected whenever (x, r) is not the 
closest neighbor of P'''(xnew, 7"now), cf. Eq. IIA.ll|l . We define the distance between 
vertices using a simple Euclidean norm ||(x, t) — (x„ew, Tnow)!! = (x — Xncw)^/-^^ + 
(r - w) 



The acceptance ratio function ljA.3|l is then 



n = . 

det A(5p) 
where Wj is given by Eq. (TO? 



detA(5;+i) 



i-U)a 



(A. 14) 



The scheme (|A.14|l in the dilute regime typically yields a an order of magnitude 
gain in efficiency, as compared to the scheme (jA.12|l . 



Appendix A. 2. 3. Shifting the worm two-point vertices Clearly, the scheme HA.14|I 
should be supplemented by an update which allows for an easy change of the closest 
neighbor. The following update performs the task: 

5p = {...,P(x,r), ...}t^{...,P(x',r'), ...} = S'^ . (A.15) 

We select x' among the nearest neighbors of the lattice site x with equal probabilities 

and r' uniformly in some interval around r. 

This update is self-balanced, and the acceptance ratio function TZ is simply 

~ 2 
det A(5p) 



n = 



detA(5;) 



(A.16) 
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